## This file writes table for replication data 

rm(list=ls())
library(ri)
load("d2ddata.rdata")
iter <- 100000
##Experiment 1, Copenhagen 

##Assign household, block, treatment, contact, and turnout

cph_data <- d2ddata[!is.na(d2ddata$itt_cph),]

cph_C <- cph_data$hh_id_all 
cph_B <- cph_data$list_cph_jod
cph_Z <- cph_data$itt_cph
cph_O <- cph_data$treated_cph
cph_Y <- cph_data$stemt

##Define treatment probabilities
probs <- genprobexact(cph_Z, clustvar = cph_C, blockvar = cph_B)

tab_cph <- 
  table(probs, cph_Y, cph_Z, cph_O)


##Experiment 2, Randers 

##Assign household, treatment, contact, and turnout

ran_data <- d2ddata[!is.na(d2ddata$itt_ran),]

ran_C <- ran_data$hh_id_all 
ran_Z <- ran_data$itt_ran
ran_O <- ran_data$treated_ran
ran_Y <- ran_data$stemt

##Define treatment probabilities
probs <- genprobexact(ran_Z, clustvar = ran_C)

tab_ran <- 
  table(probs, ran_Y, ran_Z, ran_O)

## run RI with and without clusters to find factor to adjust with

##Define treatment probabilities
probs <- genprobexact(cph_Z, clustvar = cph_C, blockvar = cph_B)

##Define approximate permutation matrix
perms_cph     <-
  genperms(cph_Z, clustvar = cph_C, blockvar = cph_B, maxiter = iter)
perms_cph_no_clust <- 
  genperms(cph_Z, blockvar = cph_B, maxiter = iter)
##Estimate ITT 

itt_cph   <- estate(cph_Y, cph_Z, prob = probs)

##gen hypothesized PO's
Ys_cph    <- genouts(cph_Y, cph_Z, ate = itt_cph)

distout   <- gendist(Ys_cph, perms_cph, prob = probs)
##Summarize distribution
dispdist(distout, ate = itt_cph)

scale_cph <-
  dispdist(distout, ate = itt_cph)$sd / 
  dispdist(gendist(Ys_cph, perms_cph_no_clust, prob = probs), ate = itt_cph)$sd 


##Define treatment probabilities
probs <- genprobexact(ran_Z, clustvar = ran_C)

##Define approximate permutation matrix
perms_ran     <-
  genperms(ran_Z, clustvar = ran_C, maxiter = iter)
perms_ran_no_clust <- 
  genperms(ran_Z, maxiter = iter)
##Estimate ITT 

itt_ran   <- estate(ran_Y, ran_Z, prob = probs)

##gen hypothesized PO's
Ys_ran    <- genouts(ran_Y, ran_Z, ate = itt_ran)

distout   <- gendist(Ys_ran, perms_ran, prob = probs)
##Summarize distribution
dispdist(distout, ate = itt_ran)

scale_ran <-
  dispdist(distout, ate = itt_ran)$sd / 
  dispdist(gendist(Ys_ran, perms_ran_no_clust, prob = probs), ate = itt_ran)$sd 

rep_tabs <- 
  list(tab_cph,
       tab_ran,
       scale_cph,
       scale_ran)

save(rep_tabs, file = "rep_tabs.rdata")